Code
# load packages 
library(tidyverse)
library(gtsummary)
library(gt)
library(tidymodels)
library(censored)
library(paletteer)
library(survey)
library(patchwork)
theme_set(theme_light())
# paletteer_d("colorBlindness::PairedColor12Steps")
col1 = "#FF9932"; col2 = "#65CCFF"
col1 = "#CC5151"; col2 = "#422CB2"
# paletteer_d("ggthemes::colorblind")
col_vec = c("#000000FF", "#009E73FF", "#CC79A7FF", "#E69F00FF", "#D55E00FF", "#56B4E9FF", "#0072B2FF")

pa_old = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_AC_C_D.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAY, ~mean(.x, na.rm = TRUE))) 

days = readRDS(here::here("data", "accelerometry", "inclusion_c_d.rds"))

wear = 
  days %>% 
  mutate(non_zero_vals = non_na_vals - zero_vals) %>% 
  filter(include_day) %>% 
  group_by(SEQN) %>% 
  summarize(wear_mins = mean(non_zero_vals)) %>% 
  mutate(SEQN = as.character(SEQN))

demo = readRDS(here::here("data", "demographics", "processed", "subset_C_D_tidy.rds")) 

NHANES 2003-2006

Survey Weighted Physical Activity Summaries by Sex

Code
df_small = pa_old %>% 
  left_join(demo) %>% 
  left_join(wear) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         wear_mins,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa, 
         mod,
         vig,
         mvpa_bout,
         vig_bout,
         starts_with("x"),
         missing_ac = na,
         median,
         logmean, 
         logmedian, 
         logtotal) 
labs = c("SEQN", "Data release cycle", "Gender",  "Wear minutes",
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "MVPA minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <5", "Time [5,10]","Time (10,25]","Time (25,50]", "Time (50,75]", "Time (75,100]", "Time (100,150]", "Time (150,200]", "Time (200,500]", 
         "Time (500,1000]", "Time (1000,2000]", "Time (2000,4000]","Time (4000,6000]","Time >6000", "Missing AC",
         "Median AC", "log10 mean AC", "log10 median AC", "log10 total AC") 

names(labs) = colnames(df_small)

df_small =
  df_small %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 8,968

1

Female
N = 4,699

1

Male
N = 4,269

1

p-value

2
Wear minutes 929.1 (123.6) 916.2 (119.5) 943.3 (126.4) <0.001
Age 45.9 (17.3) 46.4 (17.6) 45.3 (16.9) 0.022
Total AC 260,231.6 (183,008.2) 231,075.7 (200,451.9) 292,324.9 (155,469.4) <0.001
Mean AC 180.7 (127.1) 160.5 (139.2) 203.0 (108.0) <0.001
MVPA minutes 22.7 (28.3) 16.9 (29.3) 29.1 (25.8) <0.001
Moderate minutes 21.7 (22.4) 16.0 (18.8) 27.9 (24.4) <0.001
Vigorous minutes 1.0 (10.6) 0.9 (14.1) 1.2 (4.2) <0.001
MVPA (10 min bout) 4.2 (19.5) 3.8 (25.3) 4.6 (9.7) <0.001
Vig (10 min bout) 0.5 (9.6) 0.5 (12.9) 0.4 (2.9) <0.001
Time <5 895.3 (128.2) 899.8 (124.0) 890.4 (132.6) 0.028
Time [5,10] 28.7 (9.6) 29.5 (9.6) 27.7 (9.5) <0.001
Time (10,25] 51.8 (16.1) 53.2 (16.2) 50.3 (16.0) <0.001
Time (25,50] 52.3 (15.5) 53.3 (15.4) 51.1 (15.5) <0.001
Time (50,75] 36.7 (10.9) 37.3 (10.7) 35.9 (11.0) <0.001
Time (75,100] 28.8 (8.6) 29.4 (8.4) 28.1 (8.7) <0.001
Time (100,150] 44.1 (13.7) 45.1 (13.0) 42.9 (14.4) <0.001
Time (150,200] 33.2 (10.3) 34.2 (10.2) 32.1 (10.3) <0.001
Time (200,500] 113.2 (36.9) 116.5 (37.4) 109.5 (36.0) <0.001
Time (500,1000] 80.5 (35.5) 79.1 (35.1) 82.0 (35.9) 0.012
Time (1000,2000] 52.2 (33.5) 45.0 (28.3) 60.2 (36.9) <0.001
Time (2000,4000] 19.1 (18.1) 14.2 (13.4) 24.5 (20.9) <0.001
Time (4000,6000] 3.1 (9.7) 2.2 (11.7) 4.0 (6.7) <0.001
Time >6000 1.0 (9.9) 0.8 (13.0) 1.2 (4.2) <0.001
Missing AC 502.3 (133.2) 501.0 (141.3) 517.6 (3.4)
    Missing 8,959 4,691 4,268
Median AC 5.8 (88.8) 5.7 (121.7) 6.0 (16.7) 0.001
log10 mean AC 0.9 (0.2) 0.8 (0.2) 0.9 (0.2) <0.001
log10 median AC 0.2 (0.4) 0.2 (0.3) 0.3 (0.4) 0.007
log10 total AC 1,247.6 (336.2) 1,220.2 (319.1) 1,277.8 (351.6) <0.001
1

Mean (SD)

2

Design-based KruskalWallis test

Medians

Code
# unweighted, use median
df_small %>% 
  tbl_summary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight),
    statistic = list(
      all_continuous() ~ "{median} ({IQR})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 8,968

1

Female
N = 4,671

1

Male
N = 4,297

1

p-value

2
Wear minutes 920.4 (158.3) 904.0 (155.7) 939.6 (157.1) <0.001
Age 45.0 (35.0) 44.0 (34.0) 46.0 (35.0) <0.001
Total AC 225,022.3 (166,694.6) 203,857.9 (131,907.2) 259,047.7 (204,421.5) <0.001
Mean AC 156.3 (115.8) 141.6 (91.6) 179.9 (142.0) <0.001
MVPA minutes 13.9 (24.9) 9.9 (17.3) 20.8 (33.5) <0.001
Moderate minutes 13.6 (24.0) 9.7 (16.8) 20.3 (31.9) <0.001
Vigorous minutes 0.0 (0.3) 0.0 (0.1) 0.0 (0.5) <0.001
MVPA (10 min bout) 0.0 (3.3) 0.0 (1.8) 0.2 (5.0) <0.001
Vig (10 min bout) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) <0.001
Time <5 902.8 (179.0) 909.0 (171.4) 894.9 (190.3) <0.001
Time [5,10] 27.4 (12.7) 28.3 (12.9) 26.5 (12.3) <0.001
Time (10,25] 49.7 (20.9) 50.7 (21.2) 48.6 (20.7) <0.001
Time (25,50] 50.4 (19.9) 51.3 (20.1) 49.3 (20.0) <0.001
Time (50,75] 35.4 (14.1) 36.0 (14.1) 34.7 (14.2) <0.001
Time (75,100] 28.0 (11.0) 28.6 (11.0) 27.3 (11.0) <0.001
Time (100,150] 43.0 (17.2) 44.0 (17.5) 42.0 (16.9) <0.001
Time (150,200] 32.7 (13.6) 33.6 (13.6) 31.7 (13.4) <0.001
Time (200,500] 110.9 (48.4) 113.7 (49.8) 107.4 (46.3) <0.001
Time (500,1000] 74.6 (48.2) 73.3 (46.5) 76.5 (50.7) <0.001
Time (1000,2000] 42.7 (40.9) 37.5 (34.4) 49.8 (49.8) <0.001
Time (2000,4000] 12.7 (20.8) 9.4 (15.3) 18.2 (26.4) <0.001
Time (4000,6000] 0.5 (2.7) 0.3 (1.3) 1.0 (4.5) <0.001
Time >6000 0.0 (0.3) 0.0 (0.1) 0.0 (0.5) <0.001
Missing AC 529.0 (41.5) 544.0 (45.8) 519.0 (3.5) 0.3
    Missing 8,957 4,663 4,294
Median AC 0.1 (3.4) 0.0 (2.4) 0.3 (4.9) <0.001
log10 mean AC 0.8 (0.3) 0.8 (0.3) 0.9 (0.4) <0.001
log10 median AC 0.0 (0.3) 0.0 (0.3) 0.1 (0.4) <0.001
log10 total AC 1,216.0 (460.9) 1,188.1 (420.0) 1,253.0 (506.0) <0.001
1

Median (IQR)

2

Wilcoxon rank sum test

Stratified by wave

Code
df_svy =
  df_small %>% 
  filter(data_release_cycle == 3) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 4,503

1

Female
N = 2,349

1

Male
N = 2,154

1

p-value

2
Wear minutes 926 (121) 915 (118) 938 (124) <0.001
Age 46 (17) 46 (18) 45 (17) 0.069
Total AC 261,985 (222,063) 233,780 (262,098) 292,741 (162,426) <0.001
Mean AC 182 (154) 162 (182) 203 (113) <0.001
MVPA minutes 23 (33) 17 (38) 29 (26) <0.001
Moderate minutes 22 (24) 16 (22) 28 (24) <0.001
Vigorous minutes 1 (15) 1 (20) 1 (5) <0.001
MVPA (10 min bout) 5 (26) 4 (35) 5 (10) <0.001
Vig (10 min bout) 1 (13) 1 (18) 0 (3) <0.001
Time <5 896 (127) 898 (124) 894 (131) 0.4
Time [5,10] 28 (10) 30 (10) 27 (9) <0.001
Time (10,25] 51 (16) 53 (16) 49 (16) <0.001
Time (25,50] 52 (16) 53 (16) 50 (15) <0.001
Time (50,75] 36 (11) 37 (11) 35 (11) <0.001
Time (75,100] 29 (9) 30 (9) 28 (9) <0.001
Time (100,150] 44 (14) 46 (13) 43 (15) <0.001
Time (150,200] 33 (10) 35 (10) 32 (10) <0.001
Time (200,500] 114 (37) 118 (37) 110 (35) <0.001
Time (500,1000] 80 (35) 79 (35) 82 (35) 0.038
Time (1000,2000] 52 (34) 44 (28) 60 (37) <0.001
Time (2000,4000] 19 (18) 14 (14) 25 (21) <0.001
Time (4000,6000] 3 (13) 2 (16) 4 (7) <0.001
Time >6000 1 (14) 1 (18) 1 (5) <0.001
Missing AC



    95 1 (24%) 1 (28%) 0 (0%)
    512 0 (4.6%) 0 (0%) 0 (29%)
    519 0 (11%) 0 (0%) 0 (71%)
    529 1 (43%) 1 (51%) 0 (0%)
    563 1 (18%) 1 (21%) 0 (0%)
    Missing 4,500 2,346 2,154
Median AC 7 (125) 8 (172) 6 (13) 0.11
log10 mean AC 1 (0) 1 (0) 1 (0) 0.007
log10 median AC 0 (0) 0 (0) 0 (0) 0.2
log10 total AC 1,248 (335) 1,225 (321) 1,272 (348) 0.007
1

Mean (SD); n (%)

2

Design-based KruskalWallis test

Code
df_svy =
  df_small %>% 
  filter(data_release_cycle == 4) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 4,465

1

Female
N = 2,350

1

Male
N = 2,115

1

p-value

2
Wear minutes 932 (126) 917 (121) 949 (129) <0.001
Age 46 (17) 47 (18) 46 (17) 0.2
Total AC 258,490 (133,323) 228,413 (110,014) 291,907 (148,198) <0.001
Mean AC 180 (93) 159 (76) 203 (103) <0.001
MVPA minutes 22 (22) 16 (16) 29 (26) <0.001
Moderate minutes 21 (21) 16 (15) 28 (24) <0.001
Vigorous minutes 1 (3) 1 (2) 1 (4) <0.001
MVPA (10 min bout) 4 (8) 3 (7) 4 (9) <0.001
Vig (10 min bout) 0 (2) 0 (2) 0 (2) 0.038
Time <5 895 (129) 902 (124) 887 (134) 0.021
Time [5,10] 29 (10) 30 (9) 28 (10) <0.001
Time (10,25] 52 (16) 53 (16) 51 (16) 0.001
Time (25,50] 53 (15) 53 (15) 52 (16) 0.043
Time (50,75] 37 (11) 37 (11) 37 (11) 0.046
Time (75,100] 29 (9) 29 (8) 28 (9) 0.005
Time (100,150] 44 (13) 45 (13) 43 (13) 0.003
Time (150,200] 33 (10) 34 (10) 32 (10) <0.001
Time (200,500] 112 (37) 115 (37) 109 (37) <0.001
Time (500,1000] 81 (36) 79 (36) 82 (36) 0.13
Time (1000,2000] 53 (34) 46 (29) 61 (36) <0.001
Time (2000,4000] 19 (18) 14 (13) 24 (21) <0.001
Time (4000,6000] 3 (6) 2 (5) 4 (6) <0.001
Time >6000 1 (3) 1 (2) 1 (4) <0.001
Missing AC



    501 3 (46%) 3 (47%) 0 (0%)
    519 0 (3.6%) 0 (0%) 0 (100%)
    537 1 (12%) 1 (12%) 0 (0%)
    551 1 (18%) 1 (19%) 0 (0%)
    582 1 (11%) 1 (12%) 0 (0%)
    662 1 (9.5%) 1 (9.8%) 0 (0%)
    Missing 4,459 2,344 2,115
Median AC 5 (15) 4 (10) 6 (20) 0.003
log10 mean AC 1 (0) 1 (0) 1 (0) <0.001
log10 median AC 0 (0) 0 (0) 0 (0) 0.006
log10 total AC 1,248 (337) 1,215 (317) 1,284 (355) <0.001
1

Mean (SD); n (%)

2

Design-based KruskalWallis test

Survey-Weighted Physical Activity Summaries by Sex and Age

Code
survey_design =
  df_small %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
mvpa_bout
(18,30] 4 (3,4) 6 (5,7) <0.001
(30,40] 4 (3,4) 5 (4,5) 0.01738
(40,50] 4 (3,5) 5 (4,6) 0.27386
(50,60] 3 (2,4) 4 (3,5) 0.03474
(60,70] 3 (2,4) 4 (3,5) 0.38607
(70,80] 7 (-3,18) 2 (-8,3) 0.32701
vig_bout
(18,30] 0 (0,1) 0 (0,1) <0.001
(30,40] 0 (0,1) 1 (0,1) 0.01738
(40,50] 0 (0,0) 1 (1,1) 0.27386
(50,60] 0 (0,0) 0 (0,1) 0.03474
(60,70] 0 (0,0) 0 (0,0) 0.38607
(70,80] 3 (-3,8) 0 (-5,0) 0.32701
x0_5
(18,30] 914 (904,923) 877 (868,888) 0.69216
(30,40] 879 (871,887) 855 (847,868) 0.49795
(40,50] 871 (858,885) 870 (857,885) 0.03857
(50,60] 878 (868,887) 886 (876,903) 0.04307
(60,70] 916 (903,928) 934 (922,949) 0.21096
(70,80] 949 (931,966) 984 (967,997) 0.32759
x5_10
(18,30] 28 (27,28) 26 (25,27) 0.69216
(30,40] 28 (27,29) 26 (25,27) 0.49795
(40,50] 30 (29,31) 28 (27,29) 0.03857
(50,60] 31 (30,32) 29 (28,30) 0.04307
(60,70] 31 (30,32) 30 (29,30) 0.21096
(70,80] 32 (31,33) 30 (28,31) 0.32759
x10_25
(18,30] 50 (49,52) 48 (47,49) <0.001
(30,40] 52 (50,54) 49 (48,51) 0.00583
(40,50] 54 (52,56) 51 (49,53) 0.93900
(50,60] 55 (53,57) 51 (49,54) 0.31344
(60,70] 53 (52,55) 52 (50,54) 0.00617
(70,80] 56 (53,58) 52 (49,53) <0.001
x25_50
(18,30] 51 (50,52) 50 (49,51) <0.001
(30,40] 53 (51,55) 51 (50,53) 0.00583
(40,50] 54 (53,56) 52 (51,54) 0.93900
(50,60] 55 (53,56) 52 (50,54) 0.31344
(60,70] 53 (52,55) 51 (49,53) 0.00617
(70,80] 54 (52,56) 51 (49,52) <0.001
x50_75
(18,30] 36 (35,36) 35 (35,36) 0.00273
(30,40] 38 (37,39) 36 (35,38) 0.01111
(40,50] 38 (37,39) 37 (36,38) 0.00889
(50,60] 39 (38,40) 37 (36,38) 0.02781
(60,70] 37 (36,38) 35 (34,36) 0.04774
(70,80] 37 (36,39) 35 (34,36) 0.02226
x75_100
(18,30] 28 (27,28) 28 (27,29) 0.00273
(30,40] 30 (29,31) 28 (27,29) 0.01111
(40,50] 30 (29,31) 29 (28,30) 0.00889
(50,60] 31 (30,31) 29 (28,30) 0.02781
(60,70] 29 (28,30) 27 (26,28) 0.04774
(70,80] 29 (28,30) 27 (26,27) 0.02226
x100_150
(18,30] 43 (42,44) 43 (42,44) 0.00216
(30,40] 46 (44,47) 43 (42,45) 0.02779
(40,50] 46 (44,48) 44 (43,46) 0.01928
(50,60] 47 (46,48) 44 (43,45) 0.01127
(60,70] 45 (44,47) 41 (40,43) 0.07634
(70,80] 45 (44,47) 40 (38,41) 0.00789
x150_200
(18,30] 32 (31,33) 32 (31,33) 0.00216
(30,40] 35 (34,35) 33 (32,34) 0.02779
(40,50] 35 (34,36) 32 (31,34) 0.01928
(50,60] 36 (35,37) 33 (32,34) 0.01127
(60,70] 35 (34,36) 31 (30,32) 0.07634
(70,80] 34 (33,36) 30 (28,31) 0.00789
x200_500
(18,30] 109 (106,112) 110 (107,114) 0.07747
(30,40] 119 (116,121) 114 (111,118) 0.12478
(40,50] 121 (118,125) 111 (107,114) 0.10130
(50,60] 125 (121,128) 113 (110,117) 0.01986
(60,70] 121 (118,124) 107 (103,110) 0.01089
(70,80] 111 (106,117) 97 (91,101) 0.00596
x500_1e_03
(18,30] 79 (76,82) 85 (82,89) 0.07747
(30,40] 87 (83,90) 92 (88,96) 0.12478
(40,50] 88 (84,91) 86 (83,89) 0.10130
(50,60] 85 (81,88) 84 (81,88) 0.01986
(60,70] 74 (71,77) 73 (70,77) 0.01089
(70,80] 57 (53,61) 58 (54,62) 0.00596
x1e_03_2e_03
(18,30] 50 (48,52) 67 (65,71) 0.54919
(30,40] 54 (51,57) 74 (71,78) 0.16870
(40,50] 52 (50,55) 66 (64,70) 0.14382
(50,60] 45 (43,48) 59 (56,62) 0.00387
(60,70] 35 (32,37) 43 (41,47) 0.00318
(70,80] 22 (20,24) 28 (26,31) 0.00730
x2e_03_4e_03
(18,30] 18 (16,19) 31 (30,34) 0.54919
(30,40] 17 (16,18) 31 (30,34) 0.16870
(40,50] 17 (16,18) 27 (26,30) 0.14382
(50,60] 13 (11,14) 20 (19,22) 0.00387
(60,70] 9 (8,10) 14 (13,15) 0.00318
(70,80] 6 (5,7) 8 (7,9) 0.00730
x4e_03_6e_03
(18,30] 3 (2,3) 6 (5,7) 0.58355
(30,40] 3 (2,3) 5 (4,5) 0.06021
(40,50] 2 (2,3) 4 (3,5) 0.15626
(50,60] 2 (1,2) 3 (3,4) 0.00150
(60,70] 1 (1,1) 2 (2,3) <0.001
(70,80] 3 (-2,8) 1 (-4,1) <0.001
x6e_03_inf
(18,30] 1 (1,1) 2 (1,2) 0.58355
(30,40] 1 (1,1) 1 (1,2) 0.06021
(40,50] 1 (1,1) 1 (1,2) 0.15626
(50,60] 0 (0,0) 1 (1,1) 0.00150
(60,70] 0 (0,0) 0 (0,0) <0.001
(70,80] 3 (-3,8) 0 (-5,0) <0.001
total
(18,30] 249414 (240645,258183) 338721 (329952,352897) 0.87103
(30,40] 259369 (251096,267642) 344614 (336341,357723) 0.03300
(40,50] 256144 (247892,264396) 313436 (305184,328690) 0.13275
(50,60] 227431 (216765,238097) 273748 (263082,288199) 0.00158
(60,70] 188786 (179801,197772) 213610 (204624,223538) <0.001
(70,80] 182335 (109907,254763) 154885 (82457,163860) <0.001
mean
(18,30] 173 (167,179) 235 (229,245) 0.87103
(30,40] 180 (174,186) 239 (234,248) 0.03300
(40,50] 178 (172,184) 218 (212,228) 0.13275
(50,60] 158 (151,165) 190 (183,200) 0.00158
(60,70] 131 (125,137) 148 (142,155) <0.001
(70,80] 127 (76,177) 108 (57,114) <0.001
mvpa
(18,30] 21 (19,23) 38 (37,42) 0.72753
(30,40] 20 (18,21) 37 (35,39) 0.01495
(40,50] 19 (18,21) 32 (31,35) 0.01013
(50,60] 14 (13,16) 24 (22,26) <0.001
(60,70] 10 (9,12) 16 (14,17) <0.001
(70,80] 12 (1,22) 9 (-2,10) <0.001
mod
(18,30] 20 (18,22) 37 (35,40) 0.72753
(30,40] 19 (18,20) 35 (34,38) 0.01495
(40,50] 19 (17,20) 31 (29,33) 0.01013
(50,60] 14 (12,16) 23 (21,25) <0.001
(60,70] 10 (9,11) 15 (14,17) <0.001
(70,80] 9 (4,13) 8 (4,10) <0.001
vig
(18,30] 1 (1,1) 2 (1,2) 0.67430
(30,40] 1 (1,1) 1 (1,2) 0.05279
(40,50] 1 (1,1) 1 (1,2) 0.00163
(50,60] 0 (0,0) 1 (1,1) <0.001
(60,70] 0 (0,0) 0 (0,0) <0.001
(70,80] 3 (-3,9) 0 (-6,0) <0.001

Medians

Code
survey_design =
  df_small %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>% 
  filter(!is.na(age_cat))
calc_by_age =
  function(var, df) {
    df %>% 
      select(v = all_of(var), age_cat, gender) %>% 
      group_by(age_cat, gender) %>% 
      summarize(across(v, ~median(.x, na.rm = TRUE))) %>% 
      rename(median = contains("v")) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(median)) 

small %>% 
  select(metric, age_cat, Female, Male) %>% 
  group_by(metric) %>% 
  gt() %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = Female > Male 
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = Female < Male 
    )
  )
age_cat Female Male
mvpa_bout
(18,30] 0.000000e+00 2.000000e+00
(30,40] 0.000000e+00 1.000000e+00
(40,50] 0.000000e+00 5.357143e-01
(50,60] 0.000000e+00 0.000000e+00
(60,70] 0.000000e+00 0.000000e+00
(70,80] 0.000000e+00 0.000000e+00
vig_bout
(18,30] 0.000000e+00 0.000000e+00
(30,40] 0.000000e+00 0.000000e+00
(40,50] 0.000000e+00 0.000000e+00
(50,60] 0.000000e+00 0.000000e+00
(60,70] 0.000000e+00 0.000000e+00
(70,80] 0.000000e+00 0.000000e+00
x0_5
(18,30] 9.228571e+02 8.660000e+02
(30,40] 8.790000e+02 8.396000e+02
(40,50] 8.710000e+02 8.587143e+02
(50,60] 8.750000e+02 8.814167e+02
(60,70] 9.082857e+02 9.322857e+02
(70,80] 9.570000e+02 9.750000e+02
x5_10
(18,30] 2.620000e+01 2.500000e+01
(30,40] 2.680000e+01 2.500000e+01
(40,50] 2.860000e+01 2.600000e+01
(50,60] 2.900000e+01 2.750000e+01
(60,70] 2.973214e+01 2.781667e+01
(70,80] 3.166667e+01 2.824286e+01
x10_25
(18,30] 4.783333e+01 4.633333e+01
(30,40] 5.016667e+01 4.757143e+01
(40,50] 5.142857e+01 4.942857e+01
(50,60] 5.314286e+01 4.966667e+01
(60,70] 5.171429e+01 4.924286e+01
(70,80] 5.442857e+01 4.992857e+01
x25_50
(18,30] 4.820000e+01 4.814286e+01
(30,40] 5.142857e+01 5.000000e+01
(40,50] 5.320000e+01 5.080000e+01
(50,60] 5.400000e+01 5.100000e+01
(60,70] 5.214286e+01 4.824286e+01
(70,80] 5.342857e+01 4.942857e+01
x50_75
(18,30] 3.420000e+01 3.420000e+01
(30,40] 3.633333e+01 3.600000e+01
(40,50] 3.680000e+01 3.558571e+01
(50,60] 3.800000e+01 3.571429e+01
(60,70] 3.714286e+01 3.330952e+01
(70,80] 3.742857e+01 3.400000e+01
x75_100
(18,30] 2.733333e+01 2.760000e+01
(30,40] 2.885714e+01 2.814286e+01
(40,50] 2.966667e+01 2.800000e+01
(50,60] 2.983333e+01 2.816667e+01
(60,70] 2.930952e+01 2.561905e+01
(70,80] 2.933333e+01 2.630952e+01
x100_150
(18,30] 4.171429e+01 4.200000e+01
(30,40] 4.471429e+01 4.328571e+01
(40,50] 4.542857e+01 4.333333e+01
(50,60] 4.614286e+01 4.315476e+01
(60,70] 4.638095e+01 4.000000e+01
(70,80] 4.557143e+01 4.042857e+01
x150_200
(18,30] 3.120000e+01 3.200000e+01
(30,40] 3.440000e+01 3.300000e+01
(40,50] 3.471429e+01 3.253571e+01
(50,60] 3.614286e+01 3.284524e+01
(60,70] 3.550000e+01 3.057143e+01
(70,80] 3.500000e+01 3.000000e+01
x200_500
(18,30] 1.090000e+02 1.094286e+02
(30,40] 1.160000e+02 1.155000e+02
(40,50] 1.217500e+02 1.121667e+02
(50,60] 1.240000e+02 1.130000e+02
(60,70] 1.200714e+02 1.044286e+02
(70,80] 1.104286e+02 9.821429e+01
x500_1e_03
(18,30] 7.314286e+01 8.233333e+01
(30,40] 8.300000e+01 9.100000e+01
(40,50] 8.483333e+01 8.561905e+01
(50,60] 8.271429e+01 8.161905e+01
(60,70] 6.960000e+01 7.138095e+01
(70,80] 4.866667e+01 5.714286e+01
x1e_03_2e_03
(18,30] 4.250000e+01 5.880000e+01
(30,40] 4.800000e+01 6.750000e+01
(40,50] 4.614286e+01 6.021429e+01
(50,60] 4.050000e+01 5.185714e+01
(60,70] 2.771429e+01 3.775714e+01
(70,80] 1.471429e+01 2.230952e+01
x2e_03_4e_03
(18,30] 1.342857e+01 2.833333e+01
(30,40] 1.300000e+01 2.633333e+01
(40,50] 1.271429e+01 2.300000e+01
(50,60] 8.285714e+00 1.742857e+01
(60,70] 4.714286e+00 9.000000e+00
(70,80] 2.166667e+00 4.000000e+00
x4e_03_6e_03
(18,30] 7.142857e-01 3.500000e+00
(30,40] 5.714286e-01 2.666667e+00
(40,50] 4.285714e-01 1.535714e+00
(50,60] 2.857143e-01 6.904762e-01
(60,70] 1.428571e-01 2.500000e-01
(70,80] 0.000000e+00 7.142857e-02
x6e_03_inf
(18,30] 0.000000e+00 4.285714e-01
(30,40] 0.000000e+00 1.666667e-01
(40,50] 0.000000e+00 0.000000e+00
(50,60] 0.000000e+00 0.000000e+00
(60,70] 0.000000e+00 0.000000e+00
(70,80] 0.000000e+00 0.000000e+00
total
(18,30] 2.209614e+05 3.241940e+05
(30,40] 2.389592e+05 3.258449e+05
(40,50] 2.353852e+05 3.032826e+05
(50,60] 2.145245e+05 2.588716e+05
(60,70] 1.707527e+05 2.015814e+05
(70,80] 1.222527e+05 1.412056e+05
mean
(18,30] 1.534454e+02 2.251347e+02
(30,40] 1.659439e+02 2.262812e+02
(40,50] 1.634619e+02 2.106129e+02
(50,60] 1.489753e+02 1.797719e+02
(60,70] 1.185783e+02 1.399871e+02
(70,80] 8.489772e+01 9.805944e+01
mvpa
(18,30] 1.442857e+01 3.500000e+01
(30,40] 1.366667e+01 3.116667e+01
(40,50] 1.400000e+01 2.657143e+01
(50,60] 8.400000e+00 1.942857e+01
(60,70] 4.732143e+00 9.571429e+00
(70,80] 2.285714e+00 3.857143e+00
mod
(18,30] 1.400000e+01 3.316667e+01
(30,40] 1.357143e+01 3.000000e+01
(40,50] 1.360000e+01 2.564286e+01
(50,60] 8.400000e+00 1.921429e+01
(60,70] 4.714286e+00 9.500000e+00
(70,80] 2.285714e+00 3.803571e+00
vig
(18,30] 0.000000e+00 4.285714e-01
(30,40] 0.000000e+00 1.666667e-01
(40,50] 0.000000e+00 0.000000e+00
(50,60] 0.000000e+00 0.000000e+00
(60,70] 0.000000e+00 0.000000e+00
(70,80] 0.000000e+00 0.000000e+00

Stratified by wave

Code
survey_design =
  df_small %>%
  filter(data_release_cycle == 3) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 3) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
mvpa_bout
(18,30] 4 (3,6) 6 (4,7) 0.02386
(30,40] 3 (2,4) 5 (4,6) <0.001
(40,50] 4 (3,5) 5 (4,7) 0.19637
(50,60] 3 (2,4) 5 (3,6) 0.07122
(60,70] 4 (3,6) 3 (2,4) 0.34611
(70,80] 12 (-8,31) 2 (-17,3) 0.35342
vig_bout
(18,30] 1 (0,1) 1 (0,1) 0.02386
(30,40] 0 (0,1) 0 (0,1) <0.001
(40,50] 0 (0,0) 1 (1,1) 0.19637
(50,60] 0 (0,0) 0 (0,1) 0.07122
(60,70] 0 (0,0) 0 (0,0) 0.34611
(70,80] 5 (-5,16) 0 (-10,0) 0.35342
x0_5
(18,30] 916 (902,931) 877 (863,895) 0.89457
(30,40] 874 (863,885) 858 (847,878) 0.52443
(40,50] 870 (849,890) 863 (843,882) 0.08212
(50,60] 875 (866,883) 900 (891,917) 0.26054
(60,70] 910 (893,926) 953 (937,972) 0.58360
(70,80] 956 (934,977) 981 (960,1005) 0.32919
x5_10
(18,30] 27 (26,28) 26 (25,27) 0.89457
(30,40] 28 (26,30) 25 (24,26) 0.52443
(40,50] 30 (29,32) 28 (27,30) 0.08212
(50,60] 31 (29,32) 28 (26,30) 0.26054
(60,70] 31 (30,32) 30 (28,31) 0.58360
(70,80] 32 (30,34) 29 (27,30) 0.32919
x10_25
(18,30] 50 (48,51) 48 (46,50) <0.001
(30,40] 51 (49,54) 47 (45,49) 0.24128
(40,50] 54 (52,57) 52 (49,54) 0.67543
(50,60] 55 (53,58) 49 (46,52) 0.01238
(60,70] 54 (52,56) 51 (49,53) <0.001
(70,80] 55 (51,58) 50 (47,53) 0.05672
x25_50
(18,30] 51 (49,52) 49 (48,51) <0.001
(30,40] 52 (50,55) 50 (47,51) 0.24128
(40,50] 55 (52,58) 53 (50,55) 0.67543
(50,60] 56 (54,58) 49 (47,53) 0.01238
(60,70] 54 (52,55) 49 (47,50) <0.001
(70,80] 53 (50,56) 50 (47,52) 0.05672
x50_75
(18,30] 36 (34,37) 35 (34,36) 0.11573
(30,40] 37 (36,39) 35 (34,37) 0.02764
(40,50] 38 (36,40) 37 (35,39) 0.06691
(50,60] 39 (38,41) 35 (34,37) 0.04701
(60,70] 38 (36,39) 33 (32,34) 0.07384
(70,80] 36 (34,38) 34 (33,36) 0.02859
x75_100
(18,30] 28 (27,29) 28 (27,29) 0.11573
(30,40] 30 (28,31) 28 (27,29) 0.02764
(40,50] 30 (28,32) 29 (28,30) 0.06691
(50,60] 31 (30,32) 28 (27,29) 0.04701
(60,70] 30 (29,30) 26 (25,27) 0.07384
(70,80] 29 (27,30) 27 (25,28) 0.02859
x100_150
(18,30] 43 (41,44) 43 (42,45) 0.15436
(30,40] 46 (44,48) 43 (41,44) 0.03387
(40,50] 47 (44,49) 45 (43,47) 0.12130
(50,60] 48 (47,49) 42 (41,44) 0.01598
(60,70] 46 (45,48) 40 (38,41) 0.00808
(70,80] 44 (42,47) 40 (38,42) 0.03169
x150_200
(18,30] 32 (31,33) 33 (32,34) 0.15436
(30,40] 35 (34,36) 32 (31,34) 0.03387
(40,50] 36 (34,37) 33 (32,35) 0.12130
(50,60] 37 (36,38) 33 (32,34) 0.01598
(60,70] 36 (34,37) 30 (29,31) 0.00808
(70,80] 34 (32,36) 30 (28,32) 0.03169
x200_500
(18,30] 109 (104,113) 112 (108,118) 0.18842
(30,40] 122 (119,125) 113 (110,118) 0.10493
(40,50] 124 (119,129) 113 (108,118) 0.28985
(50,60] 127 (122,131) 112 (108,116) 0.00586
(60,70] 123 (118,127) 104 (99,108) <0.001
(70,80] 109 (102,116) 101 (93,106) 0.03914
x500_1e_03
(18,30] 77 (73,82) 85 (81,92) 0.18842
(30,40] 90 (85,95) 92 (87,98) 0.10493
(40,50] 87 (82,92) 87 (82,91) 0.28985
(50,60] 84 (80,87) 83 (79,87) 0.00586
(60,70] 74 (68,79) 70 (65,75) <0.001
(70,80] 56 (51,60) 60 (55,66) 0.03914
x1e_03_2e_03
(18,30] 49 (46,52) 65 (62,71) 0.30178
(30,40] 55 (52,59) 76 (73,83) 0.13340
(40,50] 50 (47,54) 66 (63,72) 0.46801
(50,60] 43 (39,46) 58 (54,61) <0.001
(60,70] 34 (29,39) 40 (35,45) <0.001
(70,80] 20 (18,22) 29 (27,33) 0.08577
x2e_03_4e_03
(18,30] 19 (17,21) 30 (27,32) 0.30178
(30,40] 17 (15,19) 34 (32,38) 0.13340
(40,50] 16 (14,17) 27 (25,31) 0.46801
(50,60] 13 (11,15) 20 (18,21) <0.001
(60,70] 10 (8,12) 13 (11,14) <0.001
(70,80] 6 (4,7) 8 (7,9) 0.08577
x4e_03_6e_03
(18,30] 3 (2,4) 6 (5,7) 0.94802
(30,40] 2 (1,3) 5 (4,6) 0.02987
(40,50] 2 (2,3) 4 (3,5) 0.42978
(50,60] 2 (1,3) 4 (3,5) <0.001
(60,70] 1 (0,2) 2 (1,2) <0.001
(70,80] 5 (-4,14) 1 (-8,2) 0.03305
x6e_03_inf
(18,30] 1 (1,2) 2 (1,2) 0.94802
(30,40] 1 (0,1) 1 (1,2) 0.02987
(40,50] 0 (0,1) 1 (1,2) 0.42978
(50,60] 0 (0,1) 1 (1,1) <0.001
(60,70] 0 (0,0) 0 (0,0) <0.001
(70,80] 5 (-5,16) 0 (-10,0) 0.03305
total
(18,30] 254386 (240503,268269) 332115 (318232,350459) 0.49035
(30,40] 262405 (249548,275261) 356634 (343778,378669) 0.01703
(40,50] 247547 (234130,260964) 313775 (300358,334544) 0.42035
(50,60] 224743 (211646,237840) 270976 (257879,286446) <0.001
(60,70] 190957 (177227,204686) 197972 (184243,212781) <0.001
(70,80] 210066 (73864,346269) 157540 (21337,172983) 0.01666
mean
(18,30] 177 (167,186) 231 (221,243) 0.49035
(30,40] 182 (173,191) 248 (239,263) 0.01703
(40,50] 172 (163,181) 218 (209,232) 0.42035
(50,60] 156 (147,165) 188 (179,199) <0.001
(60,70] 133 (123,142) 137 (128,148) <0.001
(70,80] 146 (51,240) 109 (15,120) 0.01666
mvpa
(18,30] 23 (20,27) 37 (33,41) 0.57917
(30,40] 20 (17,22) 40 (37,44) 0.01839
(40,50] 18 (16,20) 32 (30,36) 0.09202
(50,60] 15 (12,17) 23 (21,26) <0.001
(60,70] 11 (9,13) 14 (12,16) <0.001
(70,80] 16 (-4,35) 9 (-11,10) 0.02074
mod
(18,30] 22 (19,25) 35 (32,38) 0.57917
(30,40] 19 (16,21) 39 (36,43) 0.01839
(40,50] 18 (15,20) 31 (29,34) 0.09202
(50,60] 14 (12,16) 23 (20,25) <0.001
(60,70] 11 (9,13) 14 (12,16) <0.001
(70,80] 10 (2,18) 9 (0,10) 0.02074
vig
(18,30] 1 (1,2) 2 (1,2) 0.16255
(30,40] 1 (0,1) 1 (1,2) 0.02064
(40,50] 0 (0,1) 1 (1,2) 0.01938
(50,60] 0 (0,1) 1 (1,1) <0.001
(60,70] 0 (0,0) 0 (0,0) <0.001
(70,80] 6 (-5,17) 0 (-11,0) 0.04391
Code
survey_design =
  df_small %>%
  filter(data_release_cycle == 4) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 4) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
mvpa_bout
(18,30] 3 (2,4) 6 (5,7) 0.00501
(30,40] 4 (3,5) 4 (3,5) 0.74898
(40,50] 4 (3,5) 4 (3,5) 0.82775
(50,60] 3 (2,4) 4 (3,5) 0.23840
(60,70] 2 (2,3) 4 (4,6) 0.07196
(70,80] 3 (2,4) 2 (1,3) 0.38992
vig_bout
(18,30] 0 (0,1) 0 (0,1) 0.00501
(30,40] 0 (0,1) 1 (0,1) 0.74898
(40,50] 0 (0,1) 1 (1,1) 0.82775
(50,60] 0 (0,0) 0 (0,1) 0.23840
(60,70] 0 (0,0) 0 (0,0) 0.07196
(70,80] 0 (0,0) 0 (0,0) 0.38992
x0_5
(18,30] 911 (900,923) 877 (865,890) 0.52682
(30,40] 884 (874,895) 851 (841,870) 0.72103
(40,50] 873 (855,891) 877 (859,898) 0.27731
(50,60] 880 (863,898) 873 (856,899) 0.06383
(60,70] 921 (903,939) 917 (898,936) 0.24140
(70,80] 940 (912,969) 987 (958,999) 0.21847
x5_10
(18,30] 28 (27,29) 26 (25,27) 0.52682
(30,40] 29 (28,30) 28 (27,29) 0.72103
(40,50] 29 (28,31) 27 (26,29) 0.27731
(50,60] 30 (29,32) 30 (28,31) 0.06383
(60,70] 31 (29,32) 30 (28,31) 0.24140
(70,80] 32 (30,34) 31 (29,32) 0.21847
x10_25
(18,30] 51 (50,53) 48 (46,50) 0.00132
(30,40] 53 (51,55) 52 (50,53) 0.00376
(40,50] 53 (51,56) 50 (48,53) 0.78499
(50,60] 55 (53,57) 53 (51,56) 0.59360
(60,70] 53 (50,56) 53 (50,56) 0.61409
(70,80] 57 (54,59) 53 (50,56) 0.00473
x25_50
(18,30] 51 (50,52) 50 (49,51) 0.00132
(30,40] 54 (51,56) 53 (51,55) 0.00376
(40,50] 54 (52,56) 52 (50,54) 0.78499
(50,60] 54 (52,56) 54 (52,57) 0.59360
(60,70] 53 (50,56) 53 (50,56) 0.61409
(70,80] 56 (53,58) 52 (49,54) 0.00473
x50_75
(18,30] 36 (35,36) 36 (35,37) 0.01027
(30,40] 38 (37,39) 37 (36,39) 0.23841
(40,50] 38 (37,39) 36 (35,38) 0.06667
(50,60] 38 (37,39) 38 (37,39) 0.30456
(60,70] 37 (35,39) 36 (34,39) 0.27720
(70,80] 39 (36,41) 35 (33,37) 0.29542
x75_100
(18,30] 28 (27,28) 28 (27,29) 0.01027
(30,40] 30 (28,31) 29 (28,31) 0.23841
(40,50] 30 (29,31) 28 (27,30) 0.06667
(50,60] 30 (29,31) 30 (29,31) 0.30456
(60,70] 29 (27,30) 28 (26,30) 0.27720
(70,80] 30 (29,31) 27 (25,28) 0.29542
x100_150
(18,30] 43 (42,44) 42 (41,44) 0.00390
(30,40] 45 (44,47) 44 (43,47) 0.46062
(40,50] 45 (44,47) 43 (41,45) 0.08025
(50,60] 45 (44,47) 45 (43,47) 0.34405
(60,70] 44 (42,46) 43 (40,45) 0.69794
(70,80] 47 (44,49) 39 (37,41) 0.08751
x150_200
(18,30] 32 (31,33) 32 (31,33) 0.00390
(30,40] 34 (33,35) 33 (32,35) 0.46062
(40,50] 34 (33,36) 32 (30,34) 0.08025
(50,60] 35 (34,36) 33 (32,35) 0.34405
(60,70] 34 (33,36) 32 (31,34) 0.69794
(70,80] 35 (33,37) 29 (27,30) 0.08751
x200_500
(18,30] 110 (106,113) 107 (104,112) 0.25684
(30,40] 115 (111,119) 115 (111,120) 0.65398
(40,50] 118 (113,123) 109 (104,114) 0.20979
(50,60] 123 (118,127) 114 (109,119) 0.94003
(60,70] 119 (115,123) 109 (105,115) 0.86345
(70,80] 114 (104,123) 94 (84,98) 0.05127
x500_1e_03
(18,30] 81 (77,84) 84 (81,89) 0.25684
(30,40] 84 (79,88) 92 (87,97) 0.65398
(40,50] 88 (83,94) 86 (80,89) 0.20979
(50,60] 86 (81,91) 85 (80,91) 0.94003
(60,70] 75 (71,78) 76 (72,82) 0.86345
(70,80] 59 (53,66) 56 (49,61) 0.05127
x1e_03_2e_03
(18,30] 51 (48,54) 69 (66,75) 0.80340
(30,40] 53 (49,57) 72 (68,76) 0.65742
(40,50] 55 (52,57) 66 (63,70) 0.19569
(50,60] 48 (44,52) 60 (56,65) 0.89445
(60,70] 35 (32,38) 47 (43,50) 0.52245
(70,80] 24 (20,28) 28 (24,32) 0.03411
x2e_03_4e_03
(18,30] 16 (15,17) 33 (32,37) 0.80340
(30,40] 17 (15,18) 29 (27,31) 0.65742
(40,50] 18 (17,19) 28 (26,31) 0.19569
(50,60] 13 (11,15) 21 (19,24) 0.89445
(60,70] 9 (7,10) 15 (13,16) 0.52245
(70,80] 7 (5,8) 7 (6,9) 0.03411
x4e_03_6e_03
(18,30] 2 (2,3) 6 (6,8) 0.45105
(30,40] 3 (2,4) 4 (3,5) 0.52543
(40,50] 3 (2,3) 4 (3,4) 0.23567
(50,60] 2 (1,2) 3 (2,4) 0.58007
(60,70] 1 (1,1) 3 (2,4) 0.41695
(70,80] 1 (0,1) 1 (0,2) 0.00244
x6e_03_inf
(18,30] 1 (1,1) 2 (2,2) 0.45105
(30,40] 1 (1,1) 1 (1,2) 0.52543
(40,50] 1 (1,1) 1 (1,2) 0.23567
(50,60] 0 (0,0) 1 (1,1) 0.58007
(60,70] 0 (0,0) 0 (0,1) 0.41695
(70,80] 0 (0,0) 0 (0,0) 0.00244
total
(18,30] 244468 (233741,255194) 345691 (334965,367519) 0.59801
(30,40] 256310 (246017,266603) 331941 (321648,344872) 0.47110
(40,50] 264732 (254030,275434) 313103 (302401,335408) 0.16915
(50,60] 230002 (213559,246445) 276450 (260007,300438) 0.56800
(60,70] 186868 (174751,198986) 228355 (216237,240040) 0.16732
(70,80] 151700 (133850,169550) 152349 (134500,161934) <0.001
mean
(18,30] 170 (162,177) 240 (233,255) 0.59801
(30,40] 178 (171,185) 231 (223,239) 0.47110
(40,50] 184 (176,191) 217 (210,233) 0.16915
(50,60] 160 (148,171) 192 (181,209) 0.56800
(60,70] 130 (121,138) 159 (150,167) 0.16732
(70,80] 105 (93,118) 106 (93,112) <0.001
mvpa
(18,30] 19 (18,21) 40 (39,45) 0.94371
(30,40] 20 (18,22) 33 (32,36) 0.26065
(40,50] 21 (19,23) 32 (30,36) 0.04576
(50,60] 14 (12,17) 24 (22,28) 0.11235
(60,70] 10 (8,12) 17 (16,20) 0.02303
(70,80] 7 (5,9) 9 (7,11) <0.001
mod
(18,30] 18 (17,20) 38 (37,43) 0.94371
(30,40] 19 (18,21) 32 (30,34) 0.26065
(40,50] 20 (18,22) 31 (29,34) 0.04576
(50,60] 14 (12,17) 23 (21,27) 0.11235
(60,70] 10 (8,11) 17 (15,19) 0.02303
(70,80] 7 (5,9) 8 (6,10) <0.001
vig
(18,30] 1 (1,1) 2 (2,2) 0.33546
(30,40] 1 (1,1) 1 (1,2) 0.90075
(40,50] 1 (1,1) 1 (1,2) 0.02948
(50,60] 0 (0,0) 1 (1,1) 0.02907
(60,70] 0 (0,0) 0 (0,1) 0.01403
(70,80] 0 (0,0) 0 (0,0) 0.00142

Plots

Distributions

Code
pa_old %>% 
  left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(mvpa, mvpa_bout, vig, vig_bout, mean)) %>% 
  mutate(name = factor(name, labels = c("Mean AC", "MVPA (min)", "MVPA (10 min bout min)", "VPA (min)", "VPA (10 min bout min)"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2)) + 
  labs(x = "Value", y = "Density")

We see that the data are very highly skewed, so we can Winsorize at the 99th percentile:

Code
pa_old_win = 
  pa_old %>% 
    left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  mutate(across(c(mvpa, mod, vig, mean, total, mvpa_bout, vig_bout),
                ~DescTools::Winsorize(.x, quantile(.x, probs = c(0, 0.99)))))

pa_old_win %>% 
  pivot_longer(cols = c(mvpa, mvpa_bout, vig, vig_bout, mean)) %>% 
  mutate(name = factor(name, labels = c("Mean AC", "MVPA (min)", "MVPA (10 min bout min)", "VPA (min)", "VPA (10 min bout min)"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2)) + 
  labs(x = "Value", y = "Density")

Smoothed, unweighted

Code
pa_old_win %>%
  pivot_longer(cols = c(contains("bout"), vig, mvpa)) %>% 
  mutate(name = factor(name, labels = c("MVPA", "MVPA 10 min bout", "VPA", "VPA 10 min bout"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code
# pa_old_win %>%
#   pivot_longer(cols = c("mvpa", "vig", "mod")) %>% 
#   mutate(name = factor(name, labels = c("Mod", "MVPA", "VPA"))) %>% 
#   ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
#   geom_jitter(alpha = .1, size = .7, width = .25) +
#   geom_smooth(se = F) +
#   facet_grid(name~data_release_cycle, scales = "free_y") + 
#   labs(x = "Age", y = "Minutes", color = "")+
#   scale_color_manual(values = c(col1, col2)) +
#   theme(legend.position = c(0.9, 0.9),
#         legend.title = element_blank())

pa_old_win %>%
  pivot_longer(cols = starts_with("x")) %>% 
  mutate(name = factor(
    name,
    levels = c(
      "x0_5",
      "x5_10",
      "x10_25",
      "x25_50",
      "x50_75",
      "x75_100" ,
      "x100_150",
      "x150_200",
      "x200_500" ,
      "x500_1e_03",
      "x1e_03_2e_03",
      "x2e_03_4e_03",
      "x4e_03_6e_03",
      "x6e_03_inf"
    )
  )) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Smoothed, weighted

Code
survey_design =
  df_small %>%
  mutate(across(c(mvpa, mod, vig, mean, total, starts_with("x"), mvpa_bout, vig_bout),
                ~DescTools::Winsorize(.x, quantile(.x, probs = c(0, 0.99))))) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), mvpa, mod, vig, mean, total) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  mutate(metric = factor(metric, labels = c("Mean AC", "MPA (min)", "MVPA (min)", "MVPA 10 min bout (min)", "Total AC", "VPA (min)", "VPA 10 min bout (min)"))) %>%
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code
means_df =
  map_dfr(.x = pa_old %>% select(starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
   mutate(metric = factor(
    metric,
    levels = c(
      "x0_5",
      "x5_10",
      "x10_25",
      "x25_50",
      "x50_75",
      "x75_100" ,
      "x100_150",
      "x150_200",
      "x200_500" ,
      "x500_1e_03",
      "x1e_03_2e_03",
      "x2e_03_4e_03",
      "x4e_03_6e_03",
      "x6e_03_inf"
    )
  )) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Investigating outliers

Code
outliers = 
  pa_old %>% 
  left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  arrange(desc(total)) %>% 
  slice(1:10) 

# pa_old %>% 
#   arrange(desc(mvpa)) %>% 
#   slice(1:10) %>% 
#   select(SEQN, total, mean, mvpa)

# by day 
# pa_byday = readRDS(here::here("data", "accelerometry", 
#                               "summarized", "pa_summary_AC_C_D.rds")) %>% 
#   mutate(SEQN = as.character(SEQN)) 
# 
# pa_byday %>% 
#   filter(SEQN == "26729") %>% 
#   select(SEQN, PAXDAY, mod, vig, mvpa)

minute_level = readRDS(here::here("data", "accelerometry", "minute_level", "NHANES_1440_C_PAXINTEN.rds")) %>% 
  bind_rows(readRDS(here::here("data", "accelerometry", "minute_level", "NHANES_1440_D_PAXINTEN.rds")))

for(sub in outliers$SEQN){
  temp = 
    minute_level %>% 
    filter(SEQN == sub)
  
  mvpa = round(outliers %>% filter(SEQN == sub) %>% pull(mvpa), 0)
  age = demo %>% filter(SEQN == sub) %>% pull(age_in_years_at_screening)
  sex = demo %>% filter(SEQN == sub) %>% pull(gender)
  p = 
    temp %>% 
    pivot_longer(cols = starts_with("min")) %>% 
    mutate(m = sub(".*min\\_", "", name) %>% as.numeric()) %>%
    ggplot(aes(x = m, y = value)) + 
    geom_line() + 
    facet_grid(PAXDAY~.) + 
    geom_hline(aes(yintercept = 2020), color = "lightblue", linetype = "dashed") + 
    scale_x_continuous(breaks=seq(0,1440,120), labels =seq(0,24,2)) + 
    labs(x = "Hour of Day", y = "Counts", title = paste0("SEQN = ", sub, " MVPA = ", mvpa),
         subtitle = paste0("Age = ", age, " Sex = ", sex))
  
  print(p)
}

NHANES 2011-2014

Code
# load packages 


pa = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_AC.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE))) 

days = read_csv(here::here("data", "accelerometry", "inclusion_summary.csv.gz"))

wear_gh = 
  days %>% 
  group_by(SEQN) %>% 
  summarize(wear_mins = mean(non_flag_wear)) %>% 
  mutate(SEQN = as.character(SEQN))

demo_gh = readRDS(here::here("data", "demographics",  "covariates_accel_mortality_df.rds")) %>% 
  select(-contains("peak"), -contains("total_"))

Survey Weighted Physical Activity Summaries by Sex

Code
df_small_gh = pa %>% 
  left_join(demo_gh) %>% 
  left_join(wear_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         wear_mins,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x"),
         missing_ac = na) 
labs = c("SEQN", "Data release cycle", "Gender",  "Wear minutes",
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "MVPA minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <100", "Time [100,500]","Time (500,1000]","Time (1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]",  "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Missing AC") 

names(labs) = colnames(df_small_gh)

df_small_gh =
  df_small_gh %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_gh %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 9,298

1

Female
N = 4,876

1

Male
N = 4,422

1

p-value

2
Wear minutes 1,176 (137) 1,177 (135) 1,174 (140) 0.4
Age 47 (17) 48 (18) 47 (17) <0.001
Total AC 2,557,250 (825,267) 2,648,375 (817,869) 2,456,752 (821,806) <0.001
Mean AC 1,778 (574) 1,841 (569) 1,708 (572) <0.001
MVPA minutes 299 (120) 314 (118) 282 (119) <0.001
Moderate minutes 122 (41) 125 (40) 118 (43) <0.001
Vigorous minutes 137 (84) 149 (86) 125 (81) <0.001
MVPA (10 min bout) 99 (84) 108 (85) 89 (82) <0.001
Vig (10 min bout) 32 (44) 35 (45) 28 (41) <0.001
Time <100 574 (105) 580 (102) 568 (108) <0.001
Time [100,500] 129 (39) 121 (35) 137 (40) <0.001
Time (500,1000] 119 (29) 111 (26) 127 (29) <0.001
Time (1000,2000] 160 (34) 155 (31) 166 (35) <0.001
Time (2000,3000] 112 (25) 112 (24) 112 (27) 0.5
Time (3000,4000] 91 (24) 92 (23) 90 (25) 0.015
Time (4000,5000] 77 (25) 79 (24) 75 (26) <0.001
Time (5000,6000] 62 (25) 65 (25) 59 (26) <0.001
Time (6000,7000] 45 (23) 48 (23) 42 (23) <0.001
Time (7000,8000] 30 (20) 32 (20) 27 (19) <0.001
Time (8000,9000] 18 (14) 20 (15) 15 (14) <0.001
Time >9000 23 (27) 25 (28) 20 (25) <0.001
Missing AC 20 (21) 22 (22) 18 (20) 0.11
    Missing 6,809 3,736 3,073
1

Mean (SD)

2

Design-based KruskalWallis test

Stratified by wave

Code
df_svy =
  df_small_gh %>% 
  filter(data_release_cycle == 7) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 4,609

1

Female
N = 2,390

1

Male
N = 2,219

1

p-value

2
Wear minutes 1,178 (136) 1,180 (131) 1,176 (141) 0.3
Age 47 (17) 48 (17) 46 (17) <0.001
Total AC 2,565,566 (818,326) 2,655,689 (812,943) 2,468,463 (813,162) <0.001
Mean AC 1,784 (569) 1,846 (565) 1,716 (566) <0.001
MVPA minutes 299 (119) 315 (119) 282 (117) <0.001
Moderate minutes 122 (41) 126 (41) 118 (42) <0.001
Vigorous minutes 137 (83) 149 (85) 125 (79) <0.001
MVPA (10 min bout) 99 (83) 108 (85) 89 (79) <0.001
Vig (10 min bout) 32 (43) 35 (45) 27 (40) <0.001
Time <100 570 (104) 576 (101) 564 (107) 0.013
Time [100,500] 129 (39) 121 (36) 137 (40) <0.001
Time (500,1000] 119 (28) 111 (26) 128 (28) <0.001
Time (1000,2000] 162 (34) 156 (32) 168 (35) <0.001
Time (2000,3000] 113 (26) 112 (24) 113 (27) 0.7
Time (3000,4000] 91 (24) 93 (23) 90 (25) 0.071
Time (4000,5000] 77 (25) 80 (25) 75 (26) 0.003
Time (5000,6000] 62 (25) 65 (25) 59 (25) <0.001
Time (6000,7000] 45 (23) 48 (24) 42 (23) <0.001
Time (7000,8000] 30 (19) 32 (20) 27 (18) <0.001
Time (8000,9000] 18 (14) 20 (15) 15 (13) <0.001
Time >9000 23 (26) 25 (26) 20 (26) <0.001
Missing AC 19 (21) 21 (21) 18 (20) 0.3
    Missing 3,404 1,855 1,549
1

Mean (SD)

2

Design-based KruskalWallis test

Code
df_svy =
  df_small_gh %>% 
  filter(data_release_cycle == 8) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 4,689

1

Female
N = 2,487

1

Male
N = 2,202

1

p-value

2
Wear minutes 1,173 (139) 1,173 (139) 1,173 (139) >0.9
Age 47 (18) 48 (18) 47 (17) 0.015
Total AC 2,548,867 (832,210) 2,641,164 (822,799) 2,444,654 (830,654) <0.001
Mean AC 1,772 (579) 1,836 (572) 1,700 (578) <0.001
MVPA minutes 298 (121) 312 (118) 282 (122) <0.001
Moderate minutes 121 (41) 124 (39) 118 (43) <0.001
Vigorous minutes 137 (85) 149 (86) 125 (83) <0.001
MVPA (10 min bout) 99 (85) 108 (85) 90 (85) <0.001
Vig (10 min bout) 32 (45) 36 (46) 28 (43) <0.001
Time <100 578 (106) 584 (103) 572 (108) 0.002
Time [100,500] 128 (38) 120 (34) 138 (40) <0.001
Time (500,1000] 118 (29) 110 (26) 126 (30) <0.001
Time (1000,2000] 159 (33) 154 (31) 164 (34) <0.001
Time (2000,3000] 111 (25) 111 (24) 112 (26) 0.5
Time (3000,4000] 90 (24) 91 (23) 90 (25) 0.082
Time (4000,5000] 77 (25) 78 (24) 75 (26) 0.001
Time (5000,6000] 62 (25) 64 (24) 59 (26) <0.001
Time (6000,7000] 45 (23) 48 (23) 42 (24) <0.001
Time (7000,8000] 30 (20) 32 (20) 27 (19) <0.001
Time (8000,9000] 18 (15) 20 (15) 15 (14) <0.001
Time >9000 23 (27) 26 (29) 19 (25) <0.001
Missing AC 20 (21) 22 (23) 18 (20) 0.2
    Missing 3,404 1,881 1,523
1

Mean (SD)

2

Design-based KruskalWallis test

Survey-Weighted Physical Activity Summaries by Sex and Age

Code
survey_design =
  df_small_gh %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
total
(18,30] 2822818 (2775629,2870006) 2722939 (2675751,2809921) 0.01498
(30,40] 2911286 (2840893,2981680) 2690187 (2619794,2746393) <0.001
(40,50] 2814679 (2736284,2893075) 2612294 (2533899,2686557) <0.001
(50,60] 2714765 (2638819,2790712) 2437418 (2361472,2507407) <0.001
(60,70] 2465049 (2401490,2528608) 2103828 (2040269,2166496) <0.001
(70,80] 1905875 (1840890,1970860) 1706065 (1641080,1761126) <0.001
mean
(18,30] 1963 (1930,1995) 1895 (1862,1955) 0.01498
(30,40] 2024 (1975,2074) 1871 (1822,1910) <0.001
(40,50] 1957 (1902,2011) 1816 (1762,1868) <0.001
(50,60] 1887 (1834,1940) 1695 (1642,1743) <0.001
(60,70] 1713 (1669,1758) 1462 (1418,1506) <0.001
(70,80] 1324 (1279,1370) 1185 (1140,1224) <0.001
mvpa_montoye
(18,30] 332 (325,339) 315 (309,328) 0.01706
(30,40] 347 (338,356) 313 (303,321) <0.001
(40,50] 336 (325,347) 305 (294,316) <0.001
(50,60] 328 (316,340) 282 (270,292) <0.001
(60,70] 295 (286,304) 238 (229,249) <0.001
(70,80] 211 (201,221) 175 (165,184) <0.001
mod_montoye
(18,30] 121 (118,124) 120 (117,124) 0.01706
(30,40] 129 (126,132) 124 (121,127) <0.001
(40,50] 130 (127,134) 125 (122,129) <0.001
(50,60] 134 (129,139) 124 (119,128) <0.001
(60,70] 129 (125,132) 114 (110,118) <0.001
(70,80] 104 (99,108) 91 (86,94) <0.001
vig_montoye
(18,30] 173 (169,178) 158 (154,166) 0.00703
(30,40] 178 (171,185) 149 (142,156) <0.001
(40,50] 166 (158,174) 139 (131,146) <0.001
(50,60] 152 (143,160) 118 (109,125) <0.001
(60,70] 125 (118,132) 85 (78,91) <0.001
(70,80] 69 (64,74) 49 (44,54) <0.001
mvpa_montoye_bout
(18,30] 122 (118,127) 114 (109,123) 0.00703
(30,40] 130 (122,137) 107 (100,113) <0.001
(40,50] 121 (113,129) 101 (92,108) <0.001
(50,60] 114 (104,123) 84 (74,91) <0.001
(60,70] 93 (86,99) 60 (53,66) <0.001
(70,80] 49 (45,53) 33 (28,37) <0.001
vig_montoye_bout
(18,30] 45 (43,47) 40 (38,45) 0.64368
(30,40] 47 (43,51) 35 (31,38) 0.00996
(40,50] 41 (37,45) 32 (28,35) 0.05577
(50,60] 35 (30,40) 23 (18,27) 0.00356
(60,70] 25 (21,30) 13 (9,16) <0.001
(70,80] 9 (7,10) 5 (4,7) <0.001
x0_100
(18,30] 580 (573,588) 555 (547,565) 0.64368
(30,40] 557 (550,565) 541 (533,549) 0.00996
(40,50] 562 (552,573) 549 (538,558) 0.05577
(50,60] 564 (554,574) 558 (548,570) 0.00356
(60,70] 587 (578,597) 599 (589,609) <0.001
(70,80] 644 (632,656) 647 (635,658) <0.001
x100_500
(18,30] 113 (110,115) 129 (126,133) <0.001
(30,40] 113 (111,116) 132 (130,135) <0.001
(40,50] 116 (113,120) 132 (128,136) <0.001
(50,60] 116 (113,120) 137 (133,140) <0.001
(60,70] 127 (124,130) 147 (144,152) <0.001
(70,80] 147 (143,152) 162 (158,165) <0.001
x500_1e_03
(18,30] 107 (104,110) 120 (118,123) <0.001
(30,40] 108 (106,110) 125 (123,128) <0.001
(40,50] 110 (107,113) 125 (123,129) <0.001
(50,60] 110 (107,113) 130 (126,132) <0.001
(60,70] 113 (111,115) 133 (130,137) <0.001
(70,80] 121 (118,123) 134 (132,137) <0.001
x1e_03_2e_03
(18,30] 152 (149,155) 163 (160,166) 0.02708
(30,40] 154 (152,157) 167 (165,170) <0.001
(40,50] 155 (152,158) 167 (164,171) <0.001
(50,60] 157 (153,161) 169 (165,174) <0.001
(60,70] 155 (151,158) 165 (161,170) <0.001
(70,80] 158 (155,160) 164 (161,168) <0.001
x2e_03_3e_03
(18,30] 109 (107,110) 110 (109,113) 0.02708
(30,40] 111 (110,113) 113 (112,116) <0.001
(40,50] 111 (109,114) 114 (111,116) <0.001
(50,60] 114 (112,117) 115 (113,118) <0.001
(60,70] 113 (110,116) 111 (108,113) <0.001
(70,80] 112 (110,115) 112 (109,115) <0.001
x3e_03_4e_03
(18,30] 87 (86,89) 87 (85,89) 0.01478
(30,40] 92 (90,93) 91 (89,93) <0.001
(40,50] 92 (90,94) 92 (90,94) 0.00171
(50,60] 96 (94,99) 93 (90,96) <0.001
(60,70] 95 (93,98) 91 (88,93) <0.001
(70,80] 90 (87,92) 85 (82,88) <0.001
x4e_03_5e_03
(18,30] 75 (74,77) 75 (73,78) 0.01478
(30,40] 80 (79,82) 78 (76,80) <0.001
(40,50] 81 (79,83) 79 (77,81) 0.00171
(50,60] 84 (81,87) 79 (76,81) <0.001
(60,70] 82 (80,84) 74 (72,77) <0.001
(70,80] 68 (65,71) 61 (58,63) <0.001
x5e_03_6e_03
(18,30] 65 (63,66) 63 (62,66) <0.001
(30,40] 69 (67,71) 65 (63,66) 0.00584
(40,50] 69 (67,71) 64 (62,67) 0.06039
(50,60] 70 (67,73) 62 (59,65) 0.34903
(60,70] 64 (63,66) 53 (51,56) 0.08821
(70,80] 46 (44,48) 37 (35,40) 0.64263
x6e_03_7e_03
(18,30] 52 (50,53) 49 (48,51) <0.001
(30,40] 55 (53,56) 49 (47,51) 0.00584
(40,50] 53 (51,55) 47 (45,49) 0.06039
(50,60] 51 (48,54) 42 (39,44) 0.34903
(60,70] 44 (42,46) 33 (31,35) 0.08821
(70,80] 27 (25,29) 20 (18,22) 0.64263
x7e_03_8e_03
(18,30] 38 (37,39) 34 (33,36) <0.001
(30,40] 39 (37,40) 32 (31,34) <0.001
(40,50] 36 (35,38) 30 (28,32) <0.001
(50,60] 33 (31,35) 25 (23,27) <0.001
(60,70] 27 (25,29) 18 (16,19) <0.001
(70,80] 14 (13,16) 9 (8,10) <0.001
x8e_03_9e_03
(18,30] 25 (24,26) 22 (21,23) <0.001
(30,40] 25 (24,26) 20 (18,21) <0.001
(40,50] 22 (21,24) 17 (16,19) <0.001
(50,60] 19 (18,21) 14 (12,15) <0.001
(60,70] 15 (14,16) 8 (7,9) <0.001
(70,80] 7 (6,8) 4 (3,5) <0.001
x9e_03_inf
(18,30] 35 (34,37) 30 (28,32) <0.001
(30,40] 35 (31,38) 25 (22,27) <0.001
(40,50] 29 (27,32) 22 (20,25) <0.001
(50,60] 23 (21,25) 16 (14,18) <0.001
(60,70] 16 (14,18) 9 (6,10) <0.001
(70,80] 5 (5,6) 4 (3,4) <0.001

Stratified by wave

Code
survey_design =
  df_small_gh %>%
  filter(data_release_cycle == 7) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)
small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 7) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
total
(18,30] 2760827 (2703034,2818621) 2704620 (2646826,2838398) 0.38584
(30,40] 2928179 (2838858,3017500) 2738119 (2648798,2817110) 0.00115
(40,50] 2842576 (2713747,2971406) 2616878 (2488049,2736548) 0.03035
(50,60] 2794832 (2673840,2915825) 2419536 (2298543,2527517) <0.001
(60,70] 2440773 (2360832,2520714) 2144570 (2064629,2209794) <0.001
(70,80] 1890705 (1820594,1960816) 1675078 (1604967,1755619) <0.001
mean
(18,30] 1920 (1880,1960) 1882 (1842,1975) 0.38584
(30,40] 2036 (1973,2099) 1904 (1842,1960) 0.00115
(40,50] 1976 (1886,2065) 1819 (1730,1902) 0.03035
(50,60] 1943 (1859,2027) 1682 (1598,1757) <0.001
(60,70] 1696 (1640,1751) 1490 (1435,1536) <0.001
(70,80] 1314 (1265,1362) 1164 (1115,1220) <0.001
mvpa_montoye
(18,30] 323 (314,332) 310 (301,330) 0.40087
(30,40] 350 (338,363) 318 (305,329) 0.00120
(40,50] 341 (322,359) 306 (288,324) 0.03035
(50,60] 338 (319,357) 279 (260,295) <0.001
(60,70] 292 (281,304) 244 (232,253) <0.001
(70,80] 210 (199,220) 169 (158,181) <0.001
mod_montoye
(18,30] 120 (116,125) 118 (113,125) 0.40087
(30,40] 131 (127,134) 125 (121,128) 0.00120
(40,50] 132 (127,138) 126 (121,133) 0.03035
(50,60] 134 (128,141) 124 (118,130) <0.001
(60,70] 130 (126,135) 116 (111,120) <0.001
(70,80] 104 (98,110) 88 (81,93) <0.001
vig_montoye
(18,30] 165 (159,171) 155 (149,167) 0.20492
(30,40] 179 (170,188) 154 (145,162) <0.001
(40,50] 167 (155,180) 140 (127,151) 0.01984
(50,60] 162 (148,176) 114 (100,125) <0.001
(60,70] 120 (111,129) 88 (79,92) <0.001
(70,80] 67 (63,71) 46 (42,53) <0.001
mvpa_montoye_bout
(18,30] 114 (106,121) 109 (102,123) 0.20492
(30,40] 132 (123,141) 111 (101,118) <0.001
(40,50] 121 (108,134) 101 (88,112) 0.01984
(50,60] 126 (110,142) 80 (65,91) <0.001
(60,70] 87 (79,96) 61 (53,67) <0.001
(70,80] 48 (44,52) 31 (27,37) <0.001
vig_montoye_bout
(18,30] 41 (38,44) 39 (36,45) 0.59327
(30,40] 48 (43,52) 38 (33,42) 0.03672
(40,50] 40 (34,47) 31 (25,36) 0.16619
(50,60] 42 (33,50) 22 (13,26) 0.02537
(60,70] 22 (17,27) 13 (9,15) <0.001
(70,80] 8 (7,9) 5 (5,8) <0.001
x0_100
(18,30] 582 (570,593) 553 (541,570) 0.59327
(30,40] 552 (544,561) 536 (527,549) 0.03672
(40,50] 553 (539,567) 549 (535,565) 0.16619
(50,60] 559 (544,574) 557 (542,573) 0.02537
(60,70] 582 (568,596) 586 (572,603) <0.001
(70,80] 646 (630,662) 647 (631,663) <0.001
x100_500
(18,30] 114 (109,118) 130 (126,137) 0.08441
(30,40] 114 (110,117) 129 (126,133) <0.001
(40,50] 116 (110,122) 129 (123,136) 0.00704
(50,60] 114 (110,119) 138 (133,144) <0.001
(60,70] 129 (124,134) 149 (145,153) <0.001
(70,80] 149 (142,156) 162 (155,168) <0.001
x500_1e_03
(18,30] 109 (104,114) 122 (117,126) 0.08441
(30,40] 109 (106,111) 125 (122,128) <0.001
(40,50] 110 (106,114) 124 (120,128) 0.00704
(50,60] 108 (104,112) 130 (126,135) <0.001
(60,70] 115 (112,119) 135 (132,140) <0.001
(70,80] 120 (116,124) 137 (134,141) <0.001
x1e_03_2e_03
(18,30] 155 (150,160) 165 (161,169) 0.47957
(30,40] 154 (151,158) 169 (166,173) 0.00310
(40,50] 157 (154,161) 168 (164,173) 0.02780
(50,60] 156 (150,163) 171 (165,180) <0.001
(60,70] 156 (150,163) 167 (161,175) <0.001
(70,80] 156 (151,160) 168 (164,173) <0.001
x2e_03_3e_03
(18,30] 110 (108,112) 111 (109,115) 0.47957
(30,40] 112 (110,114) 114 (112,117) 0.00310
(40,50] 114 (110,117) 114 (111,119) 0.02780
(50,60] 114 (110,117) 116 (113,121) <0.001
(60,70] 115 (111,119) 111 (106,115) <0.001
(70,80] 112 (108,116) 112 (108,117) <0.001
x3e_03_4e_03
(18,30] 88 (86,90) 87 (84,91) 0.54589
(30,40] 92 (90,95) 91 (89,93) <0.001
(40,50] 94 (91,98) 93 (89,97) 0.03529
(50,60] 95 (92,99) 94 (90,98) 0.00107
(60,70] 97 (94,101) 91 (88,95) 0.00433
(70,80] 90 (86,94) 84 (79,88) 0.03401
x4e_03_5e_03
(18,30] 75 (72,78) 74 (71,78) 0.54589
(30,40] 81 (79,84) 78 (76,80) <0.001
(40,50] 83 (79,86) 79 (76,83) 0.03529
(50,60] 84 (80,88) 79 (75,83) 0.00107
(60,70] 83 (80,87) 75 (72,78) 0.00433
(70,80] 68 (64,73) 59 (55,63) 0.03401
x5e_03_6e_03
(18,30] 64 (61,66) 62 (59,66) 0.00666
(30,40] 70 (68,72) 65 (63,67) 0.01437
(40,50] 70 (66,73) 65 (61,69) 0.73753
(50,60] 71 (66,75) 61 (57,65) 0.84044
(60,70] 64 (62,66) 55 (52,57) 0.73466
(70,80] 46 (43,49) 35 (33,39) 0.89408
x6e_03_7e_03
(18,30] 50 (48,52) 48 (46,51) 0.00666
(30,40] 55 (53,58) 50 (47,52) 0.01437
(40,50] 53 (50,57) 47 (44,51) 0.73753
(50,60] 54 (49,58) 41 (37,44) 0.84044
(60,70] 43 (41,46) 33 (31,35) 0.73466
(70,80] 27 (25,28) 18 (17,21) 0.89408
x7e_03_8e_03
(18,30] 36 (35,38) 33 (32,36) 0.00119
(30,40] 39 (37,41) 33 (31,35) <0.001
(40,50] 37 (34,40) 30 (27,33) 0.00828
(50,60] 36 (32,39) 24 (21,27) <0.001
(60,70] 26 (24,28) 18 (15,19) <0.001
(70,80] 14 (13,15) 9 (8,10) 0.00982
x8e_03_9e_03
(18,30] 23 (22,25) 21 (20,23) 0.00119
(30,40] 25 (24,26) 20 (19,22) <0.001
(40,50] 23 (20,25) 18 (15,19) 0.00828
(50,60] 21 (19,24) 13 (11,15) <0.001
(60,70] 14 (12,16) 8 (7,9) <0.001
(70,80] 7 (6,7) 4 (3,5) 0.00982
x9e_03_inf
(18,30] 33 (30,35) 30 (28,33) <0.001
(30,40] 34 (30,38) 27 (24,30) <0.001
(40,50] 29 (26,32) 22 (18,25) <0.001
(50,60] 26 (22,29) 15 (11,18) <0.001
(60,70] 14 (11,17) 10 (7,11) <0.001
(70,80] 5 (4,5) 4 (3,5) <0.001
Code
survey_design =
  df_small_gh %>%
  filter(data_release_cycle == 8) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 8) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
total
(18,30] 2882509 (2808573,2956446) 2742097 (2668160,2853457) 0.00689
(30,40] 2895061 (2787659,3002462) 2639481 (2532080,2719143) 0.00645
(40,50] 2786316 (2696997,2875635) 2607407 (2518088,2692178) 0.00323
(50,60] 2631658 (2553190,2710126) 2457263 (2378795,2543668) 0.00888
(60,70] 2489119 (2393179,2585059) 2065497 (1969557,2170579) <0.001
(70,80] 1920711 (1811994,2029429) 1734750 (1626032,1810594) 0.00703
mean
(18,30] 2004 (1953,2056) 1909 (1857,1986) 0.00689
(30,40] 2013 (1939,2088) 1835 (1761,1891) 0.00645
(40,50] 1938 (1876,2000) 1813 (1751,1872) 0.00323
(50,60] 1829 (1774,1883) 1709 (1654,1768) 0.00888
(60,70] 1730 (1663,1797) 1435 (1368,1508) <0.001
(70,80] 1335 (1259,1411) 1206 (1130,1258) 0.00703
mvpa_montoye
(18,30] 341 (331,350) 321 (311,336) 0.00761
(30,40] 344 (330,358) 308 (294,319) 0.00630
(40,50] 332 (319,344) 303 (291,316) 0.00311
(50,60] 317 (303,331) 286 (272,299) 0.00916
(60,70] 298 (285,311) 233 (220,252) <0.001
(70,80] 212 (195,229) 182 (165,194) 0.00714
mod_montoye
(18,30] 122 (119,125) 122 (119,127) 0.00761
(30,40] 128 (124,132) 124 (120,128) 0.00630
(40,50] 128 (124,132) 124 (121,129) 0.00311
(50,60] 133 (127,140) 124 (118,130) 0.00916
(60,70] 127 (123,131) 111 (107,119) <0.001
(70,80] 103 (97,110) 93 (87,99) 0.00714
vig_montoye
(18,30] 181 (175,188) 162 (155,173) 0.00555
(30,40] 176 (165,188) 145 (133,154) 0.00343
(40,50] 165 (156,174) 139 (130,148) <0.001
(50,60] 141 (134,149) 122 (115,130) 0.00541
(60,70] 130 (120,141) 83 (73,94) <0.001
(70,80] 71 (61,80) 52 (42,58) 0.00887
mvpa_montoye_bout
(18,30] 131 (125,137) 118 (113,130) 0.00555
(30,40] 127 (117,138) 102 (92,112) 0.00343
(40,50] 121 (111,131) 100 (91,111) <0.001
(50,60] 101 (93,109) 88 (79,96) 0.00541
(60,70] 98 (88,108) 59 (48,69) <0.001
(70,80] 50 (42,57) 34 (27,40) 0.00887
vig_montoye_bout
(18,30] 50 (46,53) 42 (38,48) 0.93807
(30,40] 47 (40,53) 32 (26,37) 0.11891
(40,50] 41 (36,47) 33 (27,38) 0.15608
(50,60] 28 (25,31) 25 (22,31) 0.06455
(60,70] 29 (23,35) 13 (7,17) 0.00153
(70,80] 10 (7,12) 5 (3,7) 0.03433
x0_100
(18,30] 579 (569,589) 556 (547,568) 0.93807
(30,40] 562 (550,574) 546 (534,556) 0.11891
(40,50] 572 (556,588) 549 (533,558) 0.15608
(50,60] 570 (557,582) 559 (546,576) 0.06455
(60,70] 593 (580,607) 611 (598,623) 0.00153
(70,80] 642 (625,659) 647 (630,662) 0.03433
x100_500
(18,30] 112 (110,115) 127 (124,132) <0.001
(30,40] 112 (109,116) 135 (131,139) 0.00217
(40,50] 117 (113,120) 134 (131,139) <0.001
(50,60] 119 (113,124) 136 (131,140) 0.00410
(60,70] 125 (121,129) 145 (141,153) <0.001
(70,80] 146 (141,150) 161 (157,164) 0.00346
x500_1e_03
(18,30] 105 (104,106) 119 (118,122) <0.001
(30,40] 108 (105,111) 126 (123,129) 0.00217
(40,50] 110 (107,113) 127 (123,131) <0.001
(50,60] 111 (107,116) 129 (124,132) 0.00410
(60,70] 111 (108,113) 130 (128,136) <0.001
(70,80] 121 (118,124) 132 (128,137) 0.00346
x1e_03_2e_03
(18,30] 149 (148,151) 160 (159,164) 0.01716
(30,40] 154 (151,158) 165 (162,168) 0.00848
(40,50] 153 (148,158) 166 (162,171) 0.00181
(50,60] 158 (153,163) 167 (162,170) 0.02001
(60,70] 153 (149,157) 162 (158,168) <0.001
(70,80] 159 (156,163) 161 (157,165) 0.00329
x2e_03_3e_03
(18,30] 107 (105,109) 109 (107,112) 0.01716
(30,40] 111 (108,114) 113 (110,116) 0.00848
(40,50] 109 (106,112) 113 (110,115) 0.00181
(50,60] 115 (112,118) 114 (111,117) 0.02001
(60,70] 111 (107,115) 111 (107,114) <0.001
(70,80] 113 (109,116) 111 (108,115) 0.00329
x3e_03_4e_03
(18,30] 87 (85,88) 87 (85,90) 0.00961
(30,40] 91 (88,94) 91 (88,93) 0.00471
(40,50] 90 (88,92) 91 (89,94) 0.02262
(50,60] 97 (93,101) 92 (88,96) 0.42677
(60,70] 94 (90,97) 90 (87,93) <0.001
(70,80] 89 (85,93) 87 (82,90) 0.00839
x4e_03_5e_03
(18,30] 76 (74,77) 76 (74,79) 0.00961
(30,40] 80 (77,82) 78 (75,80) 0.00471
(40,50] 80 (78,82) 78 (76,81) 0.02262
(50,60] 84 (80,88) 78 (74,82) 0.42677
(60,70] 80 (77,83) 73 (70,78) <0.001
(70,80] 68 (64,72) 63 (59,66) 0.00839
x5e_03_6e_03
(18,30] 66 (64,68) 65 (63,67) 0.00350
(30,40] 68 (66,70) 64 (62,66) 0.10196
(40,50] 68 (65,70) 64 (61,66) 0.01199
(50,60] 68 (64,72) 63 (59,66) 0.22575
(60,70] 64 (62,67) 52 (49,56) 0.01183
(70,80] 46 (42,50) 39 (35,42) 0.60374
x6e_03_7e_03
(18,30] 54 (52,55) 51 (49,54) 0.00350
(30,40] 54 (51,56) 48 (46,51) 0.10196
(40,50] 52 (50,55) 46 (44,49) 0.01199
(50,60] 49 (46,52) 43 (40,45) 0.22575
(60,70] 45 (42,47) 32 (29,36) 0.01183
(70,80] 27 (24,30) 21 (18,23) 0.60374
x7e_03_8e_03
(18,30] 40 (38,41) 36 (34,38) <0.001
(30,40] 38 (36,40) 32 (30,34) <0.001
(40,50] 36 (34,38) 29 (27,32) <0.001
(50,60] 31 (29,33) 26 (24,28) <0.001
(60,70] 28 (26,31) 17 (15,20) <0.001
(70,80] 15 (13,17) 10 (8,11) <0.001
x8e_03_9e_03
(18,30] 26 (25,28) 22 (21,24) <0.001
(30,40] 25 (23,27) 19 (17,21) <0.001
(40,50] 22 (21,24) 17 (16,19) <0.001
(50,60] 18 (16,19) 14 (13,16) <0.001
(60,70] 16 (14,18) 8 (6,10) <0.001
(70,80] 7 (6,9) 4 (3,5) <0.001
x9e_03_inf
(18,30] 38 (35,40) 30 (27,33) <0.001
(30,40] 35 (30,40) 23 (18,26) <0.001
(40,50] 30 (27,33) 23 (20,27) <0.001
(50,60] 20 (18,21) 16 (15,19) <0.001
(60,70] 18 (14,22) 8 (4,9) <0.001
(70,80] 6 (5,8) 4 (2,4) 0.00206

Plots

Distributions

Code
pa %>% 
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(contains("montoye"), mean)) %>% 
  mutate(name = factor(name, labels = c("LPA", "Mean AC", "MPA", "MVPA", "MVPA 10 min bout", "Sedentary", "VPA", "VPA 10 min bout"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2))

Smoothed, unweighted

Code
pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = contains("bout")) %>% 
  mutate(name = factor(name, labels = c("MVPA 10 min Bout", "VPA 10 min Bout"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") +
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code
pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = contains("montoye") & !contains("bout")) %>% 
  mutate(name = factor(name, labels = c("Light", "Moderate", "MVPA", "Sedentary", "Vigorous"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code
pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = starts_with("x")) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Smoothed, weighted

Code
survey_design =
  df_small_gh %>%
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  filter(!grepl("x", metric)) %>% 
  # mutate(metric = factor(metric, labels = c("MVPA", "VPA"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code
results %>% 
  filter(grepl("x", metric)) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

MIMS

Code
pa_mims = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_PAXMTSM.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE))) 
Code
df_small_gh_mims = pa_mims %>% 
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         starts_with("x"),
         missing_mims = na) 
labs = c("SEQN", "Data release cycle", "Gender",  
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total MIMS", "Mean MIMS",
         "Time <2", "Time [2,5]","Time (5,10]","Time (10,15]", "Time (15,20]", "Time (20,30]", "Time (30,40]", "Time (40,50]", "Time (50,60]",  "Time >60", "Missing MIMS") 

names(labs) = colnames(df_small_gh_mims)

df_small_gh_mims =
  df_small_gh_mims %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_gh_mims %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 9,298

1

Female
N = 4,876

1

Male
N = 4,422

1

p-value

2
Age 47 (17) 48 (18) 47 (17) <0.001
Total MIMS 13,496 (3,819) 13,820 (3,766) 13,139 (3,845) <0.001
Mean MIMS 9 (3) 10 (3) 9 (3) <0.001
Time <2 592 (110) 596 (107) 587 (112) 0.003
Time [2,5] 151 (43) 142 (40) 160 (46) <0.001
Time (5,10] 185 (42) 177 (39) 195 (44) <0.001
Time (10,15] 134 (31) 132 (29) 136 (33) 0.002
Time (15,20] 110 (29) 110 (28) 109 (30) 0.065
Time (20,30] 167 (60) 173 (59) 161 (61) <0.001
Time (30,40] 74 (47) 80 (48) 67 (45) <0.001
Time (40,50] 19 (20) 21 (21) 17 (19) <0.001
Time (50,60] 4 (7) 4 (7) 4 (7) <0.001
Time >60 2 (6) 2 (5) 3 (8) <0.001
Missing MIMS 20 (21) 22 (22) 18 (20) 0.10
    Missing 6,805 3,734 3,072
1

Mean (SD)

2

Design-based KruskalWallis test

Code
survey_design =
  df_small_gh_mims %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_mims %>% select(starts_with("x"), total, mean) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh_mims %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
x0_2
(18,30] 594 (586,601) 569 (561,579) <0.001
(30,40] 571 (563,579) 555 (547,563) 0.01062
(40,50] 578 (567,588) 567 (556,577) 0.13192
(50,60] 581 (571,591) 579 (569,592) 0.77927
(60,70] 607 (597,618) 623 (613,635) 0.02304
(70,80] 667 (655,678) 678 (666,689) 0.09110
x2_5
(18,30] 134 (131,137) 152 (149,157) <0.001
(30,40] 135 (132,138) 156 (153,160) 0.01062
(40,50] 138 (135,142) 154 (151,160) 0.13192
(50,60] 138 (134,142) 161 (157,165) 0.77927
(60,70] 148 (145,152) 170 (166,175) 0.02304
(70,80] 167 (163,172) 181 (177,185) 0.09110
x5_10
(18,30] 173 (169,176) 190 (186,193) <0.001
(30,40] 174 (172,177) 196 (194,200) <0.001
(40,50] 177 (173,180) 194 (191,199) <0.001
(50,60] 178 (173,183) 199 (194,205) <0.001
(60,70] 177 (173,181) 198 (194,204) <0.001
(70,80] 183 (179,186) 192 (189,196) <0.001
x10_15
(18,30] 129 (127,131) 133 (131,136) <0.001
(30,40] 133 (130,135) 138 (136,141) <0.001
(40,50] 132 (130,135) 139 (136,143) <0.001
(50,60] 135 (132,138) 139 (136,143) <0.001
(60,70] 132 (129,136) 133 (129,136) <0.001
(70,80] 134 (131,137) 134 (132,138) <0.001
x15_20
(18,30] 105 (103,106) 105 (103,108) <0.001
(30,40] 110 (108,112) 110 (108,112) <0.001
(40,50] 111 (108,113) 111 (109,114) <0.001
(50,60] 115 (112,118) 111 (108,115) <0.001
(60,70] 114 (111,117) 109 (106,113) <0.001
(70,80] 111 (108,115) 107 (103,110) 0.00141
x20_30
(18,30] 170 (166,174) 166 (162,172) <0.001
(30,40] 182 (178,187) 172 (167,175) <0.001
(40,50] 182 (177,187) 171 (166,177) <0.001
(50,60] 186 (179,192) 168 (161,174) <0.001
(60,70] 176 (171,180) 151 (147,158) <0.001
(70,80] 136 (130,143) 118 (112,123) 0.00141
x30_40
(18,30] 95 (93,98) 87 (85,92) 0.02278
(30,40] 96 (92,99) 81 (78,85) <0.001
(40,50] 89 (85,93) 75 (71,79) 0.00277
(50,60] 82 (77,86) 64 (59,68) 0.11019
(60,70] 67 (62,71) 45 (40,48) 0.79112
(70,80] 35 (32,38) 24 (22,27) 0.83662
x40_50
(18,30] 29 (28,31) 26 (24,27) 0.02278
(30,40] 28 (26,30) 22 (19,23) <0.001
(40,50] 24 (22,26) 19 (17,21) 0.00277
(50,60] 20 (18,21) 13 (11,14) 0.11019
(60,70] 14 (12,16) 7 (6,8) 0.79112
(70,80] 5 (4,6) 3 (3,4) 0.83662
x50_60
(18,30] 6 (6,7) 6 (6,7) 0.89984
(30,40] 6 (6,7) 5 (4,6) 0.84073
(40,50] 5 (5,6) 4 (4,5) 0.75824
(50,60] 4 (3,4) 3 (2,3) 0.14261
(60,70] 2 (2,3) 2 (1,2) 0.03213
(70,80] 1 (0,1) 1 (0,1) 0.05500
x60_inf
(18,30] 3 (3,4) 4 (4,5) 0.89984
(30,40] 3 (2,3) 3 (3,4) 0.84073
(40,50] 3 (2,3) 3 (3,3) 0.75824
(50,60] 1 (1,2) 2 (2,3) 0.14261
(60,70] 1 (1,1) 1 (1,1) 0.03213
(70,80] 0 (0,0) 0 (0,0) 0.05500
total
(18,30] 14641 (14419,14863) 14423 (14202,14824) 0.21164
(30,40] 15034 (14712,15357) 14272 (13949,14528) <0.001
(40,50] 14568 (14216,14920) 13865 (13513,14219) 0.00306
(50,60] 14102 (13758,14446) 13017 (12673,13341) <0.001
(60,70] 12967 (12665,13269) 11428 (11127,11722) <0.001
(70,80] 10433 (10131,10736) 9607 (9305,9873) <0.001
mean
(18,30] 10 (10,10) 10 (10,10) 0.21164
(30,40] 10 (10,11) 10 (10,10) <0.001
(40,50] 10 (10,10) 10 (9,10) 0.00306
(50,60] 10 (10,10) 9 (9,9) <0.001
(60,70] 9 (9,9) 8 (8,8) <0.001
(70,80] 7 (7,7) 7 (6,7) <0.001
Code
pa_mims %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(total, mean))  %>%
  mutate(name = factor(name, labels = c("Mean MIMS", "Total MIMS"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Value", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Code
pa_mims %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = starts_with("x")) %>% 
  mutate(name = factor(name, levels = c("x0_2", "x2_5", "x5_10", "x10_15", "x15_20", "x20_30", "x30_40", "x40_50", "x50_60", "x60_inf"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Code
survey_design =
  df_small_gh_mims %>%
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }


means_df =
  map_dfr(.x = pa_mims %>% select(starts_with("x"), total, mean) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>% 
  filter(!grepl("x", metric)) %>% 
  mutate(metric = factor(metric, labels = c("Mean MIMS", "Total MIMS"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Value",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code
results %>% 
  filter(grepl("x", metric)) %>% 
  mutate(metric = factor(metric, levels = c("x0_2", "x2_5", "x5_10", "x10_15", "x15_20", "x20_30", "x30_40", "x40_50", "x50_60", "x60_inf"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

NYFFS

Tables

Code
demo_y = readRDS(here::here("data", "demographics", "subset_Y_tidy.rds")) 
pa_y = readRDS(here::here("data", "accelerometry", "summarized", 
                          "pa_summary_AC_paxy.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE)))

df_small_y = pa_y %>% 
  left_join(demo_y) %>% 
  select(SEQN,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_interview_weight,
         full_sample_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         contains("montoye"),
         starts_with("x")) 
labs = c("SEQN", "Gender",  
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "Sedentary minutes", "Light minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <100", "Time [100,500]","Time (500,1000]","Time (1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]",  "Time (7000,8000]", "Time (8000,9000]", "Time >9000") 

names(labs) = colnames(df_small_y)

df_small_y =
  df_small_y %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_y %>% 
  select(SEQN,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_interview_weight,
         full_sample_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa_montoye_bout,
         vig_montoye_bout,
         contains("range")) %>% 
  mutate(WTMEC4YR_norm = full_sample_mec_exam_weight/mean(full_sample_mec_exam_weight, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_interview_weight, full_sample_mec_exam_weight, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic

Overall
N = 1,370

1

Female
N = 674

1

Male
N = 696

1

p-value

2
Age 9 (4) 9 (4) 9 (4) 0.2
Total AC 3,538,697 (805,142) 3,546,168 (755,892) 3,531,471 (850,561) 0.7
Mean AC 2,460 (559) 2,465 (525) 2,455 (591) 0.7
MVPA (10 min bout) 181 (85) 184 (83) 177 (88) 0.10
Vig (10 min bout) 80 (49) 79 (46) 81 (52) 0.7
1

Mean (SD)

2

Design-based KruskalWallis test

Code
survey_design =
  df_small_y %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(3, 6, 9, 12, 15))) %>%
  mutate(weight = full_sample_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_y %>% select(contains("bout"), contains("x"), contains("montoye")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_y %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(3, 6, 9, 12, 15))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_mec_exam_weight,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2)))  %>%
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )
age_cat Female Male pval
mvpa_montoye_bout
(3,6] 243 (227,258) 233 (218,252) 0.16257
(6,9] 213 (193,232) 193 (174,205) 0.03998
(9,12] 159 (148,170) 153 (141,162) 0.30361
(12,15] 115 (100,130) 108 (93,118) 0.41696
vig_montoye_bout
(3,6] 103 (93,113) 107 (97,117) 0.16257
(6,9] 96 (83,109) 92 (80,102) 0.03998
(9,12] 68 (61,74) 70 (64,75) 0.30361
(12,15] 47 (39,54) 49 (42,55) 0.41696
x0_100
(3,6] 563 (555,571) 551 (543,560) 0.42823
(6,9] 543 (533,553) 533 (523,540) 0.58236
(9,12] 551 (544,559) 539 (531,548) 0.49074
(12,15] 578 (557,599) 581 (561,600) 0.55775
x100_500
(3,6] 67 (64,69) 72 (70,76) 0.42823
(6,9] 73 (69,77) 84 (80,88) 0.58236
(9,12] 91 (85,97) 106 (101,112) 0.49074
(12,15] 114 (109,118) 130 (125,135) 0.55775
x500_1e_03
(3,6] 65 (63,68) 72 (70,75) 0.02874
(6,9] 76 (73,79) 86 (82,89) 0.13839
(9,12] 93 (88,97) 104 (99,107) 0.06680
(12,15] 110 (107,113) 121 (118,125) 0.76776
x1e_03_2e_03
(3,6] 113 (109,118) 123 (118,127) 0.02874
(6,9] 132 (126,138) 142 (136,147) 0.13839
(9,12] 152 (146,158) 157 (152,160) 0.06680
(12,15] 160 (156,164) 165 (161,170) 0.76776
x2e_03_3e_03
(3,6] 108 (103,112) 110 (105,115) 0.01209
(6,9] 116 (110,121) 119 (113,122) 0.00154
(9,12] 118 (115,121) 116 (113,119) <0.001
(12,15] 111 (108,115) 106 (103,110) <0.001
x3e_03_4e_03
(3,6] 102 (98,105) 102 (98,105) 0.01209
(6,9] 101 (97,106) 101 (96,104) 0.00154
(9,12] 96 (94,97) 92 (90,95) <0.001
(12,15] 86 (83,89) 79 (76,82) <0.001
x4e_03_5e_03
(3,6] 94 (91,96) 91 (89,93) 0.00191
(6,9] 88 (85,91) 85 (82,88) 0.00159
(9,12] 80 (77,83) 76 (73,78) 0.00121
(12,15] 70 (67,73) 63 (61,66) <0.001
x5e_03_6e_03
(3,6] 81 (79,83) 78 (75,80) 0.00191
(6,9] 75 (73,77) 69 (67,71) 0.00159
(9,12] 66 (63,68) 62 (60,65) 0.00121
(12,15] 58 (54,61) 52 (49,55) <0.001
x6e_03_7e_03
(3,6] 66 (65,68) 63 (61,65) <0.001
(6,9] 61 (58,63) 55 (52,57) 0.01115
(9,12] 53 (51,56) 50 (47,52) 0.07327
(12,15] 46 (43,49) 42 (39,44) 0.06065
x7e_03_8e_03
(3,6] 51 (49,53) 47 (45,50) <0.001
(6,9] 47 (44,49) 42 (39,43) 0.01115
(9,12] 41 (39,43) 37 (35,39) 0.07327
(12,15] 34 (32,37) 30 (27,32) 0.06065
x8e_03_9e_03
(3,6] 37 (35,39) 35 (33,37) 0.27116
(6,9] 35 (32,37) 30 (28,31) 0.19712
(9,12] 30 (29,32) 27 (25,29) 0.44564
(12,15] 24 (22,27) 20 (18,22) 0.04788
x9e_03_inf
(3,6] 91 (84,99) 97 (90,104) 0.27116
(6,9] 93 (84,101) 94 (85,102) 0.19712
(9,12] 68 (62,74) 71 (65,76) 0.44564
(12,15] 47 (41,53) 48 (41,53) 0.04788
sed_montoye
(3,6] 901 (890,913) 913 (901,928) 0.95553
(6,9] 925 (909,941) 948 (932,958) 0.70090
(9,12] 990 (973,1006) 1008 (992,1024) 0.03792
(12,15] 1059 (1037,1082) 1091 (1069,1109) 0.00855
light_montoye
(3,6] 111 (107,115) 111 (107,115) 0.95553
(6,9] 111 (106,116) 111 (106,115) 0.70090
(9,12] 105 (103,108) 102 (100,105) 0.03792
(12,15] 95 (92,98) 87 (84,91) 0.00855
mod_montoye
(3,6] 151 (147,155) 146 (142,149) 0.10790
(6,9] 141 (137,145) 134 (130,139) 0.11076
(9,12] 127 (122,131) 120 (116,125) 0.02906
(12,15] 111 (106,116) 101 (96,105) 0.00849
vig_montoye
(3,6] 276 (264,288) 270 (258,284) 0.10790
(6,9] 262 (246,278) 246 (230,255) 0.11076
(9,12] 217 (205,229) 208 (196,217) 0.02906
(12,15] 172 (157,187) 159 (144,170) 0.00849
mvpa_montoye
(3,6] 472 (459,484) 461 (448,476) 0.04812
(6,9] 447 (430,463) 424 (408,434) <0.001
(9,12] 385 (369,401) 368 (352,382) 0.03768
(12,15] 320 (300,341) 293 (273,309) 0.04173

Smoothed, unweighted

Code
pa_y %>%
  left_join(demo_y) %>% 
  pivot_longer(cols = contains("bout")) %>% 
  mutate(name = factor(name, labels = c("MVPA", "VPA"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes in MVPA", color = "")+
  scale_y_continuous(breaks=seq(0,600,60)) + 
  scale_color_manual(values = c(col1, col2)) + 
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank()) + 
  scale_x_continuous(breaks=seq(3,15,2))

Smoothed, weighted

Code
survey_design =
  df_small_y %>%
  mutate(weight = full_sample_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_y %>% select(contains("bout")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  mutate(metric = factor(metric, labels = c("MVPA", "VPA"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_grid(. ~ metric) +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(3, 15, 2)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code
means_df =
  map_dfr(.x = pa_y %>% select(starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
    scale_x_continuous(breaks=seq(3, 15, 2)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))